Data collapse in the critical region using finite-size scaling with subleading corrections 
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We propose a treatment of the subleading corrections to finite-size scaling that preserves the 
notion of data collapse. This approach is used to extend and improve the usual Binder cumulant 
analysis. As a demonstration, we present results for the two- and three-dimensional classical Ising 
models and the two-dimensional, double-layer quantum antiferromagnet. 



True phase transitions occur only in systems with an 
infinite number of degrees of freedom. Nevertheless, it 
is possible to study critical phenomena by simulating fi- 
nite systems of increasing size and extrapolating their 
behaviour to the thermodynamic limit. The standard 
analysis relies on Fisher's finite-size scaling (FSS) hy- 
pothesis P, , which supposes that near a critical point 
every thermodynamic property of the system scales as a 
universal function of the ratio L/£, where L is the linear 
size of the system and £ is the bulk correlation length. 

The basic ansatz is that any property with critical ex- 
ponent k obeys the relation A(T,L) = L K /"f A (L/£(t)), 
which is valid for large L and small reduced temperature 
t = (T — T c )/T c . Here, T c is the critical temperature of 
the transition, and v is the critical exponent of the cor- 
relation length £ ~ \t\~ v . The function /a (a;) depends 
only on the boundary conditions and the system geome- 
try. Equivalently, one can write 

A(T,L)L-«''' = g A (tL 1 ' v ). (1) 

Equation (JTJ can be formally motivated 0, 0, by 
arguing that the diverging correlation length leaves the 
system scale invariant As t — > 0, fluctuations of the 
dynamical variables begin to look the same on all length 
scales (much greater than the lattice spacing). Conse- 
quently, any change in the system size L can be compen- 
sated by a corresponding change in temperature t. The 
two otherwise independent quantities become locked to- 
gether in a single composite variable x — tL 1 ' v . 

Graphically, this effect is undersood in terms of so- 
called data collapse 0,0. A plot of Eq. □ in reduced 
variables traces out the universal function g A (x). Such 
a plot can be useful as a data analysis tool. Given some 
set of measurements {(Tj, Li,Ai)}, the critical behaviour 
of the system can be determined by constructing the co- 
ordinates Xi — L,y u (Ti — T c )/T c and t/, = AiL i K ^ u and 
choosing the values of v, k, and T c such that the graph 
points {(xi,yi)} have the least amount of scatter 

The accuracy of such a fit depends not only on the 
quantity and quality of the data; it also depends sensi- 
tively on whether the data are sufficiently deep in the 
scaling limit (L — > oo at fixed x). What constitutes 
"deep enough," however, is model dependent and diffi- 
cult to determine. For high-accuracy studies, it is best 
not to put complete trust in Eq. ft}. Since true scaling is 



realized only asymptotically, data from all but the largest 
systems can rarely be made to fall onto a single curve. 
Worse, the data can sometimes be shoehorned into col- 
lapse with incorrect parameters. For lattice sizes that 
are accessible to computer simulation, there are almost 
always significant, non-universal corrections to FSS. To 
simply ignore them can lead to subtle systematic errors, 
unreliable fits, and overly-optimistic error estimates. 

In this Letter, we show how to incorporate sublead- 
ing corrections into the minimal-scatter optimization de- 
scribed above. We argue that the correct approach is to 
make small, size-dependent modifications to the prefac- 
tor and argument of g A {x) in Eq. ft}: 

A(T, L)L-I V = N{L)g A {tlM v - e(L)). (2) 

In other words, the features of A are both renormalized 
and shifted. In the thermodynamic limit, Af(L) — ► 1 
and e(L) — > 0; we determine the precise asymptotic form 
of these functions by renormalization group arguments. 
The key point is that y = A(T, U)L~ K I V M(L)~ X now 
behaves as a universal function of x = tL x l v — e(L). 

In the special case where g A (x) is bounded and mono- 
tonic, A(T, L)L~ K / l> becomes increasingly step-function- 
like as L — > oo. Thus, two such curves, plotted versus 
temperature for different linear sizes L and L' , will inter- 
sect (except in unusual circumstances) at a unique point 
in the vicinity of the critical temperature T c . Binder's 
spin-distribution cumulants have this useful prop- 
erty (as does, e.g., the spin stiffness ^3). ^ n the Binder 
cumulant case, k = happens to be known a priori. We 
solve formally for the point of intersection T{(L,L') un- 
der the assumption of Eq. J5J|-like scaling and determine 
its location as a function of L and n = L'/L. We prove 
that a sequence of intersection points (n fixed, L — > oo) 
converges to T c faster than L~ x l v . Fitting the n depen- 
dence of a data set offers another way of estimating T c 
and v. 

Corrections to finite-size scaling — FSS can be put on 
a (somewhat) rigourous basis by invoking the renormal- 
ization group (RG) |3| |j, L| • From the RG point of view, 
the system is characterized by a set of scaling fields {CO- 
The singular part of the free energy behaves according to 

/ s (Ci, C2, Ca, • • •) = L- d H(iL v \ C 2 ^ 2 , C 3 L^ ,...), (3) 
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where J 7 (xi,X2,X3, . . .) is a regular function of its argu- 
ments. The yi are eigenvalues of the RG transformation 
and characterize the flow of the fields under the rescaling 
L i-> L/b with b > 1; i.e., Q >-> C*& w - 

For concreteness, suppose that there is only one rele- 
vant scaling field £1 = t with a scaling exponent z/i = 
1/v > 0. Then / s admits a series expansion in the re- 
maining fields (with Q = fixed): 



-T(tL 1 /» 1 x 2 ,x 3 ,...) 



dx r - 



=o 
(4) 

ft follows that any thermodynamic quantity generated 
from Eq. wm have the form 

A(T, L) = L K l v g A {tL x l v ) + L^-^l u p A {tL l ' v ) + • • • (5) 

where = min(|j/i|). The functions gA(x) and p A (x) 
are partial derivatives of T . Their asymptotic behaviour, 
g A (x) ~ M _K and Pa(^) ~ as x — > ±oo, leads 

to the thermodynamic limit 



A(T) = lim A(T, L) 

L — > oc 



|*r(l + (const) 1*1* 



(6) 



The functions g A (x) and Pa(^) are well-behaved and 
admit series expansions g A (x) = go + g\x + \g 2 x 2 + ■ ■ ■ , 
etc. Identifying Eqs. J5J and (JSJ up to 0(x 3 ), we find 
that 7V(i) = 1 + CiL-M" and e(L) = C 2 L-^ U , where 



Ci 



g2Po - aw\ 
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and C2 = 



3? • 



(7) 



This analysis is flawed, however, in that scaling holds 
only asymptotically. The "equalities" in Eqs. (0) and (@J 
are subject to corrections analytic in L _1 . This means 

that, strictly speaking, Af(L) = (1 + G X L-^I V H )(1 + 

aiL^ 1 + a 2 L~ 2 + •••). Over some range of L values, this 
can be represented by N(L) = 1 + cL~ u , where uj has 
some effective value close to that of the true L — > 00 
exponent, min(0/z/, 1). Thus, 

A(T, L) = L K I V {\ + cL-^gAitL 1 ^ - dL~^ v ) (8) 

is the appropriate generalization of Eq. Most im- 

portant, this scaling ansatz represents an improvement 
over Eq. (J5J in that data collapse can still be engineered 
by plotting AL,- K I V {\ + cL^)" 1 versus tL l ' v - dL'*/ 1 ', 
with the addition of ui, 4>, c, and d as fitting parameters. 

Intersection point analysis — In the vicinity of a mag- 
netic phase transition, expectation values of powers of 
the magnetization obey the scaling relation (\M\ p )l = 
L-pP/v g^ p (tL 1 / v ) for integer p > 0. At large x, the scal- 
ing function g\M\ P (x) exhibits power-law behaviour with 
exponents p(3 as x — > —00 and pf3 < p(3 as x — > 00, which 
ensures that the magnetization has the correct form in 
the thermodynamic limit: (iMp^ = (-£]ff 0(-t). One 
can define a family of quotient functions 

(M 2p )l 



Qp(T,L) 



(My L 



m 

Af{L)q p {tL 1 / v - e(L)) 



(9) 



— let us call them Binder ratios — constructed in such a 
way that the leading-order L dependence outside the scal- 
ing function factors out. The temperature dependence 
appears only in the argument of the function q p (x) = 
gM2 P (x)/[gM2(x)] p . The subleading corrections enter as 
shown on the right-hand side of Eq. © ■ The two lowest 
order Binder cumulants [Eqs. (11) and (12) of Ref. H(i| 
correspond to the linear combinations U = 1 — \Qv and 
V = 1 - \Q 2 + JLQ 3 . 

In the thermodynamic limit, the Binder ratios jump 
discontinuously between Q P (T < T c ) = q p (— 00) = 1 and 
Q P (T > T c ) = q p (co). The value of g p (oo) 7^ 1 is simply 
a prefactor of the gaussian distribution of thermally ran- 
domize spins; the value Q P (T C ) — q p (0) is a non-trivial, 
universal constant. The Binder cumulants are also step 
functions, taking the values U — 2/3 and V — 8/15 in the 
ordered phase and U = V = in the disordered phase. 

We now consider the point of intersection Ti(L, L') be- 
tween two Binder ratio curves for system sizes L =/= L' . To 
start, let us assume that T{(L,nL), for some fixed ratio 
n = L' /L, converges to T c as a function of L faster than 
L~ x l v . In that case, ti{L,nL)L x ' v is a small quantity. 
Accordingly, 



Q(T i ,L)=Af(L) 



qo + qi (tiL 



(10) 



[We have dropped the p label; the subscripts here refer to 
the expansion coefficients of q(x) about x = 0.] Hence, 
the crossing criterion Q(T U L) = Q(T[, L') implies that 



U = 



-qo [Af(L) - N{L')} + gi [Af(L)e(L) - N(L')e{L')} 



qi[N{L)V-f" -N{L>){L'y/»] 



or, in the notation of Eq. (JHJ), 



(11) 



cqofl-n "\ r _ 1/w _ w 
1 - rr^ v 



( - J~ v - \ L -^y». (12) 



[In retrospect, the assumption that Ti(L, nV) — > T c faster 
than L~ x l v was justified.] Note that although Eq. |Q 
was derived with the Binder ratio in mind, it applies 
equally to any quantity A(T, L)L~ K / U obeying Eq. © 
whose scaling function is bounded and monotonic. 

Numerical results — We have applied these results 
to three numerical test cases: the two- (2D) and 
three-dimensional (3D) nearest-neighbour Ising models, 
and the double-layer quantum Heisenberg antiferromag- 
net 0. Monte Carlo data were generated for the Ising 
models using the Swendsen-Wang cluster update algo- 
rithm |l.'.i| and for the quantum antiferromagnet using 
stochastic series expansion [l^ . 

The Ising models exhibit classical, thermally-driven 
phase transitions at T Ci 2D = 2/ log(l + \/2) w 2.2691853 
and T C .3D ~ 4.511. The antiferromagnetic bilayer, on the 
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FIG. 1: (color online) The second Binder ratio Q2(T,L), computed via Monte Carlo for the 3D Ising model on lattices of size 
L = 4, 5, . . . , 14. In the left-hand panel, the solid (green) lines are guides to the eye, connecting data points computed for equal 
system sizes; the curves grow steeper with increasing L. The oval inset presents a magnified view in which the individual data 
points and their error bars (blue) are more clearly visible. The scaling plot in the right-hand panel shows the complete collapse 
of the data onto a single curve. This best fit gives T c = 4.5114(2), v = 0.625(1), and q 2 (0) = 1.60(1) [or U* = 1 - |g 2 (0) = 
0.467(3)]. The two inset panels are purely leading-order scaling plots (c = d = 0) based on (top-left) the best fit values just 
quoted and (bottom-right) the values T c = 4.500(2), v = 0.60(1) obtained by fitting Q2 versus tL 1 ^" . 



other hand, exhibits a zero-temperature, quantum phase 
transition as the interlayer coupling stren gth a is tuned 
through its critical value a c w 2.5218 [ill Il2j. For each 
of these cases, we computed the second Binder ratio over 
a fine, uniform mesh of temperature (coupling) values in 
a fixed interval containing T c (a c ). (In practice, when 
rough estimates of T c and v are available in advance, 
it makes more sense to take measurements in the range 
|x| = I tlX 1 '" < 1.) The simulations were performed for a 
range of relatively small lattice sizes. 

The data were fit to Eq. JSJ (with k = 0) using the 
the Levenberg-Marquardt nonlinear optimization algo- 
rithm [lij ]. which searched the space of parameters T c 
(a c ), v, u>, (f>, c, and d to produce the best collapse of 
the data. Uncertainties in those values were computed 
by bootstrapping the regression over the original 
Monte Carlo data. Comparison with the 2D Ising model, 
where the solution is known |l7j . serves as an important 
proof of concept: the fitted values T c — 2.26917(2) and 
v = 0.99(1) with L = 4, 5, . . . , 32 agree within statistical 
uncertainties with the exact results. The same procedure 
applied to the quantum bilayer with L = 8, 10, 12, . . . , 42 
yields a c — 2.52181(4), the most accurate value to date, 
and v = 0.715(2), which is consistent with the 3D classi- 
cal Heisenberg universality class, as expected [lij . 

The raw data for the 3D Ising model are shown in the 
left-hand panel of Fig. ^ In the right-hand panel, the 



same data are rescaled according to Eq. JHJ and collapse 
convincingly onto a single curve. The two inset panels on 
the right illustrate the failure of naive FSS: plotted in the 
conventional reduced coordinates (c = d = Q), the data 
are clearly distinguishable as a series of separate curves 
corresponding to different lattice sizes. 

The fitted value of the critical temperature T c = 
4.5114(2) is consistent with the Rosengren conjecture 
and with other reliable estimates; cf. Ref. [2QJ and Table 
19 in Ref. |^. Note that our result is nearly as accurate 
as those of Pawley et al. and Garcia et al. (viz., T r . = 
4.5115(1) in Ref. and T c = 4.51152(12) in Ref. 
which are based on simulations up to sizes L = 64 and 
L = 115, respectively. And while our value of v = 
0.625(1) is small with respect to some estimates |2ll E2I 
Eij , it appears to be in very good agreement with more re- 
cent Monte Carlo (L = 90, 100, 115) [H and Monte Carlo 
Renormalization Group (L — 64, 128, 256) |25L l26| calcu- 
lations. As an additional check, we can read off the value 
92(0) = 1.60(1) (the "y- intercept" in the right-hand panel 
of Fig. [IJ, which agrees with the value (72(0) = 1.604(1) 
computed by independent methods |2lj . 

Figure [5] plots the Binder ratio intersection points of 
the 2D Ising model and the quantum bilayer for several 
system size ratios n = L'/L. In each case, the com- 
plete set of crossing points was computed by interpolat- 
ing smooth curves between the measured Q2 values (as 
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ing limit has an interesting implication: it may be more 
profitable to collect a large quantity of easily-obtainable 
intermediate- L data [23] (so as to maximize the fit statis- 
tics) than to make a herculean effort to obtain data for 
the very largest possible lattice sizes. 



FIG. 2: (color online) Binder ratio crossing points Ti(L,L') 
plotted with respect to L (< L') for the (left) 2D Ising model 
and (right) double-layer quantum antiferromagnet. The solid 
lines depict the analytic result given in Eq. 11211 for lines of 
constant n = L'/L. From top (red) to bottom (blue), the 
ratios shown are (left) n = 8/7, n = 3/2, 2, . . . , 11/2, 6 (by 
increments of 1/2) and (right) n = 3/2, 2, . . . , 9/2, 5. 



in the left panel of Fig. 2] where the grid lines form a 
set of intersection points). The data were successfully fit 
to the scaling form given by Eqs. Hll|) and (|12l) . The 
resulting estimates T c = 2.2692(8), v = 0.99(2) and 
a c = 2.5218(2), v — 0.715(8) are somewhat less accu- 
rate than but consistent with the values determined by 
the data collapse analysis. 

Conclusions — Finite-size scaling describes how critical 
behaviour emerges from finite systems in the limit of in- 
creasingly large system size. As L — > oo, quantities in the 
critical region become universal modulo a size-dependent 
rescaling of the variables, a property which, in principle, 
can be exploited to measure unknown critical parameters 
by means of data collapse. In practice, such an analysis 
is likely to give misleading results. Still, we have shown 
that data collapse can be made a viable, high-accuracy 
analysis tool, so long as subleading corrections to finite- 
size scaling are properly taken into account. We have 
introduced a new way of expressing those corrections, 
in which two kinds of deformation — the renormalization 
cja — * NgA and the shift gA (x) — > gA(x — e) — serve as 
the fundamental deviations from leading-order FSS. The 
greater expressiveness of the scaling form typically leads 
to larger but more meaningful and reliable statistical er- 
rors on the critical parameters. 

Large simulations devour computing resources. With 
Monte Carlo, CPU time scales as L d+Z , where z is a 
characteristic dynamical exponent. For the Ising model 
in d = 3, Swendsen-Wang updates [13J give z ~ 0.75. 
This implies that it takes as much time to compute the 
single system size L = 16 as it does to compute all of 
L = 4, 5, . . . , 12. (Even worse, for models where special- 
ized update schemes are not available, z ~ 2.) The suc- 
cess of Eq. (JSJ as a fitting form well away from the seal- 
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